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ABSTRACT 

We carry out several isolated galaxy evolution simulations in a fixed dark matter 
halo gravitational potential using the new version of our N-body/Smoothed Particle 
Hydrodynamics (SPH) code GCD+. The new code allows us to more accurately model 
and follow the evolution of the gas and stellar components of the system including 
powerful supernovae feedback and its effects on the inter-stellar medium (ISM). Here 
we present the results of six simulations of an M33-sized late- type disc galaxy each with 
varying values for our model parameters which include the star formation efficiency 
(C*), the energy released per supernovae explosion (i?sN) and the energy released per 
unit time from stellar winds (-Esw)- We carry out both a pixel- by-pixel and radial 
ring analysis method for each of our galaxies comparing our results to the observed 
Schmidt-Kennicutt Law and vertical gas velocity dispersion versus radius relation 
amongst others. We find that our models with higher feedback more closely resemble 
the observations and that feedback plays a pivotal role in obtaining both the observed 
Schmidt-Kennicutt and gas velocity dispersion relations. 

Key words: galaxies: formation — galaxies: evolution — galaxies: ISM — galaxies: 
structure — galaxies: kinematics and dynamics — ISM: bubbles 
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1 INTRODUCTION 

The process of galaxy formation and evolution remains one 
of the most outstanding problems in modern astrophysics. 
From a cosmological point of view, great progress has been 
made including the development of a theory; ACDM cos- 
mology which accurately explains the present day large scale 
distribution of galaxies from the primordial small scale den- 
sity fluctuations in the early universe. The actual processes 
which govern individual galaxy formation and evolution are 
far less well understood however. One of these impo rtant 
processes is that of star formation itself. [Schmidtl (|l959l l sug- 
gested that the star formation rate (SFR) density and the 
interstellar gas density in our Galaxy were related by a sim- 
ple power law relation of the form psf r oc (pgas)", where n 
is the power law index. ISchmidd (ligsgp originally suggested 
a value for n of 2, so that the rate of star formation varies 
with the square of the gas de nsity. Stud ying the form of the 
observed star formation law, iKennicutd (19891 . Il998bl ) anal- 
ysed the SFR surface density and the gas surface density for 
a much larger sample of gala^xies and over a wider range of 



E-mail: ara2@mssl.ucl.ac.uk 



gas densities, d etermining the pow er law index n to be equal 
to 1.4 (see also lKennicuttlll998al ) . We now refer to the em- 
pirical relation between SFR and gas surface density as the 
Schmidt-Kennicutt (SK) relation. Recentl y, the SK relation 
has been studied in much more depth (e.g. Bigiel at al.ll20oj : 
Leroy et al. 2008) including observing the relation for indi- 
vidual components of the gas (HI and molecular) and as a 
function of environment. Several observat ional studies foun d 
evidence for a two -slope SK relation (e.g. lBigiel et al.ir2008l '). 
iBigiel et al.l (|2008l ) showed that a break exists at gas surface 
densities of 9 M0yr~^, below which the slope of the SFR 
steepens. This break corresponds to the transition between 
the HI and molecular phase of the gas. Theoretical mod- 
els also predict such br eaks in the SK relation. For example 
iKrumholz et al.l (|2009l ') made a theoretical model consist- 
ing of three prescriptions: molecular self-shielding, internal 
feedback and turbulent regulation in molecular clouds which 
address the physics responsible for determining the SFR on 
local scales and successfully managed to reproduce the ob- 
served trends in the SK relation. 

The SK relation is of great importance as it is used for 
the calibration of numerical codes of galaxy formation and 
evolution, with the SFR often explicitly required to obey a 
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SK type law (e.g. ISpringel fc HernquistI l2003l : iBooth et all 
l2007l : IShave fc Dalla Vecchia 2008). Some more recent stud- 
ies however, do not force the SFR to follow a SK-like rela- 
tion, and s o its adherence to this law is a prediction of the 
model (e.g. iRobertson fc Kravtsovlliooil : ISaitoh et al1l2008l : 
[Monaco et al.ll201ll '). 

Recently IXamburro et all (|2009l ) showed that the gas 
velocity dispersion is also an important factor to consider 
for any galaxy evolution model. Using HI line width maps 
obtained for 11 disc galaxies of Th e HI Nearby Galaxy Sur- 
vey (THINGS. [Walter et al1l2008h they determined that all 
of their sample galaxies showed a systematic radial decline 
in their HI velocity dispersion, with a characteristic value of 
~ 10 km s~^ by r25 — roughly the radial extent of the star 
forming disc. The radial decline in velocity dispersion was 
well correlated with the radial decline of SFR and therefore 
they suggested that high gas velocity dispersion is induced 
by supernovae feedback. 

The implications of such studies are that the gas veloc- 
ity dispersion should now be used as an additional constraint 
for any n ew numerical simul ation codes of galaxy evolution. 
Recently. iDobbs et al.l (|201ll ) performed isolated disc galaxy 
evolution simulations to investigate how the properties of the 
ISM, molecular clouds and SFR depend on the level of stellar 
feedback, comparing their results against both the observed 
SK relation and gas velocity dispersion. They found that 
supernovae were vital in reproducing the observed velocity 
dispersions and scale heights of the individu al components 
of the ISM. Recently. iPilkington et al.1 (120111 ) analysed both 
the SK and HI velocity dispersion relations for a fully consis- 
tent cosmologically simulated late-type dwarf spiral galaxy. 

In this paper we initiate a project to create an accu- 
rate and self-consistent numerical model for late-type spiral 
galaxies. To this end we h a,ve been developing our n umer- 
ical simulation code GCD+ (|Kawata fc GibsorJ l2003al ) . The 
new version of the code includes new features and improve- 
ments to a number of aspects of the modelling of galaxy 
formation and evolution, including new star formation and 
feedback recipes. The overall effect of these changes is that 
the code can now capture and accurately follow the effects 
of powerful supernovae explosions on the evolution of the 
galaxy. The details of the new model are described in the 
next section. In this study, we set up an isolated disc galaxy 
consisting of gas and stellar discs with no bulge component, 
in a static dark matter halo gravitational potential. To cal- 
ibrate our new model, we here focus on M33-/iA;e bulgeless 
galaxies. M33 is a suitably chosen late-type bulgeless spiral 
galajcy. Its close proximity in the Local Group means it has 
been a favourite choice for detailed o bservational campaigns 
(e.g. iHever et al.ll2004 IVerlev et al.ir2010 ). Its smaller size 
compared to the Milky Way also means that we can sim- 
ulate it to much higher resolution using the same number 
of parti cles. M33 also lacks a supermassive black hole at its 
centre jMerritt et al.l l200ll ) , a factor which helps simplify 
the modelling and analysis. Therefore, M33-/iA;e galaxies are 
ideal target galaxies for testing and calibrating the numerical 
modelling of the physical processes in galaxy evolution. Note 
that our aim is not to reproduce completely the observed 
properties of M33. It is likely that M33 had a recent inter- 
actio n with M31 iMcConnachie et al.l |200^ : iPutman et al.l 
I2OO9I ) and so it may not be a typical gala:xy. We focus more 



on the general observed properties for M33-sizeci disc galax- 
ies, such as the SK relation and the gas velocity dispersion. 

The paper is organised as follows. In Section 2, we de- 
scribe the galaxy model together with some new features of 
our numerical simulation code. In Section 3, we present our 
results and analyses. We summarise our findings and present 
our conclusions and planned future work in Section 4. 



2 THE CODE AND MODEL 
2.1 Nev^r version of GCD+ 

To simulate our galaxy, we use an updated ver- 
sion of our origi nal galactic chem o dynam ical evolu- 
tion code GCD+ l|Kawata fc GibsonI 12003^ ). GCD+ is 



a parallel three-dimension al N-body/SPH (|LucvI Il977l : 
iGingold fc Monaghan|[T977l ') code which can be used in stud- 
ies of galaxy formation and evolution in both a cosmological 
and isolated setting. GCD+ incorporates self-gravity, hydro- 
dynamics, radiative cooling, star formation, supernova feed- 
back, and metal enrichment. The new version of the code 
includes metal diffusion suggested by Greif et al. (2009). 
Our recent update and performance in various test suites 
will be presented in a series of our forthcoming papers 
(Kawata et al. in preparation). Here we briefly describe the 
major points of our recent upgrade. We have implemented 
a mod ern scheme of SPH suggested by iRosswog fc Price 
1120071) . We follow their artificial viscosity switch (Morris 
[19971) and artificial thermal conductivity to re solve the 
Kelvi n- Helmholtz instabi l ity ([Kawata et al.ll200^ V Follow- 
ing iSEHngeT^^^emguisl ([2002! ') ■ instead of the energy equa- 
tion, we integrate the entropy equation. As suggested by 
[Saitoh fc Makind i 200a), we have added the individual time 
step limiter, which is crucial for correctly resolving the ex- 
pansion bubbles induced by supernovae feedback (see also 
[Merlin et al] [2OI0I : [Purier fc Dalla Vecchia[[201ll'). W e also 
implement the FAST scheme ([Saitoh fc Makind[2010l ') which 
allows the use of different timesteps for integrating hydro- 
dynamics and gravity. The code now also includes adap- 
tive s oftening for both the stars and gas ([Price fc Monaghad 
[-'nul l although in this current study we only implement it 
for the gas particles (see below). 

Radiative coo ling and heating is calculated with 
CLOUDY (vOS.OO: [Perland et al.|[l99i '). We tabulate cool- 
ing and heating rates and the mean molecular weight as 
a function of redshift, metallicity, density and temp e rature 
adopting the 2005 version of the [Haardt fc Madaul ([l996[ 'l 
UV backgroun d radiation. We also add a th ermal energy 
floor following [Robertson fc Kravtsovl ([2008h to keep the 
Jeans mass higher than 2A'nb'^^p, where A^nb is the num- 
ber of neighbour particles and mp is the particle mass, to 
avoid numerical instability (see also [Bate fc Burkcrt 199'^). 
In the resolution of simulations presented in this paper (see 
below), the gas particle meets this condition around riH = 1 
cm""^. Therefore, we allow the gas particle to become a star 
particle if the density becomes greater than nn = 1 cm^'^, 
i.e. the star formation threshold density, and if their velocity 
fie ld is convergent, following the Schmidt law as described 
in [Kawata fc Gibsod ([2003a[ ): 



dp* 
~dt 



dps 
dt 



(1) 
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Table 1. Simulation parameters 



Cnfw 




^vir 








fld.ga. ^d.star 




(Mo) 


(kpc) 


(Mo) 


(M©) 


(kpc) 


(kpc) 


10 


5.0 X 10" 


163 


4.0 X 10^ 


2.0 X 10'' 


1.4 


4.0 0.3 



where C* is a dimensionless SFR parameter and tg is the 
dynamical time. 

We assume that stars are distributed according 
to the ISalpeteil l|l955h initial mass function (IMF). 
GCD+ takes into account chem ical enrichment by both 
Type II supernovae (SNe II, IWoosle y fc Weaverl 
and Type la superno vae (SNe la, llwamoto et al.l 



1995) 



199C 



iKobavashi et all bOOOH and mass loss from intermediate- 
mass stars ( van den Hoek fc GroenewegenI [l997l 'l . and fol- 
lows the chemical enrichment history of both the stellar and 
gas components of the system. The new version of GCD+ uses 
a different scheme for star formation and feedback. We now 
keep the mass of the ba ryon (gas and star) particl es com- 
pletely the same, unlike iKawata fc GibsonI (|2003b[ ) or the 
majority of SPH simulations which includ e star formation . 
Alt hough the basic strategy is sim ilar to iLia et aP ()2002l ) 
and iMartmez-Serrano et al.l (|2008l ). a slightly different im- 
plementation is adopted. Details of the methodology will 
be described in our forthcoming paper. Here we describe it 
briefly. First, every star particle formed in the simulation is 
randomly assigned a mass group ID ranging from 1 to 61 
(although 61 is chosen arbitrarily, it is a compromised se- 
lection to sample the stellar mass range and be similar to 
our resolution, i.e. number of neighbour particles). We can 
calculate that about 13 per cent of the mass is ejected by 
SNe II from a star cluster following the assumed IMF and 
IWooslev fc Weaveil (|l995l ). Therefore we set ID=l-8 parti- 
cles to be 'SNe II particles'. Each SNe II particle receives 1/8 
of the energy and metals produced by SNe II depending on 
the age and initial metallicity. After the age reaches a time 
when 8/61 of the mass is ejected (about the lifetime of a 8 
M0 star, depending on metallicity), the whole star particle 
is changed back to a gas particle. We assume each super- 
nova produces thermal energy Esn (ergs). We also assume 
that stellar winds from massive stars (> 30 Mq) produce 
thermal energy -Bsw (ergs s"^) and add this to the SNe 
particles. We explore the effect of these parameters apply- 
ing Esn = 10f^_and lO^ergs and Esw = 10^*^ and 10^'^ 
ergs s (e.g. lGibsonlll994l ). Consequently SNe II particles 
have higher temperature and become metal rich. We do not 
take into account radiative cooling until the particle turns 
back to a gas particle. However, we calculate thermal en- 
ergy for the SNe II particles following the SPH scheme and 
add the pressure from the SNe II particles to their neigh- 
bour particles, while the dynamics of the SNe II particle s 
are only affected by gravity (see also lPelupessv et"al]|2004 ). 
The metals produced by SNe II are also distributed from 
the SNe II particl es through the metal diffusion scheme of 
iGreif et al.l (|2009l ). i.e. the metal diffusion applied to the 
SNe II particles is the same as the other normal gas par- 
ticles. We apply a similar algorithm to star particles with 
different IDs. From the IMF about 30 per cent of the mass 
will be ejected after a cosmic time. Therefore star particles 
with IDs ranging from 9 to 19 will turn back into a gas 



particle depending on the age of the star particle, which 
models the mass loss from intermediate mass stars and SNe 
la feedback. For these particles, we turn on radiative cool- 
ing during their mass loss and SNe la. For example, ID 9 
particles with solar metallicity represent mass-loss from 7 
M0 to 5.6 Mq stars. When their age becomes older than 
the lifetime of a 7 Mq star, they become a 'feedback parti- 
cle'. The particle inherits the original metal abundance and 
receiv es additional metals that stars in this mass range pro- 
duce (|van den Hoek fc GroenewegenI Il997l ). which are then 
diffused via the metal diffusion scheme applied. Their ther- 
mal energy is calculated with the SPH scheme, and the ad- 
ditional pressure from these feedback particles are applied 
to their neighbouring particles (especially for the higher ID 
particles with SNe la feedback). Once the particle becomes 
older than the lifetime of a 5.6 Mq star, the particle becomes 
a normal gas particle. Due to this algorithm, the particle 
mass of stars and gas are always constant, and the mass 
resolution is kept constant. 

The main free parameters of the new version of GCD+ 
include the star formation efficiency, C, , energy per super- 
nova, -BsN, and stellar wind energy per massive star, -Esw- 
Especially these parameters control the effect of feedback 
on star formation, which is the most unknown and possibly 
most important process in galaxy formation and evolution. 
This paper explores the effect of these parameters and cal- 
ibrates these parameters against the observed properties, 
such as the SK relation and the gas velocity dispersion. 



2.2 Model galaxy setup 

We are interested in investigating the evolution of a late- type 
disc galaxy of the order of size of M33. We therefore set up 
an isolated disc galaxy which consists of gas and stellar discs 
with no bulge component in a static dark matter halo poten- 
tial. We use the standard Na varro- Frenk- White ( NFW) dark 
matter halo density profile (|Navarro et al.|[l997l ). assuming 
a standard cold dark matter (ACDM) cosmological model 
with cosmological parameters of Q,o = 0.266, flb = 0.044 
and Ho = 71kms"^Mpc"^: 



Pd^ - (1 + zo) ^^^^ ^^^^ ^ , 
where 



_ J'200 _ r 

C : ^ ) 

rs '■200 



(2) 



(3) 



and 



.2oo^l.63xlO-(^)^[|a^]-^(l+.o)-;.-kpc,(4) 

where pc is the characteristic density of the profile, r is 
the distance from the centre of the halo and is the 
scale radius. The halo mass is set to be M200 ~ 5.0 x 
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Run 1 Run 2 




Figure 1. Face on (top row) and edge on (bottom row) view of the stellar (left column) and gas (right column) distribution for our six 
simulation runs at t = 1 Gyr when we carry out our analyses. The large holes in the gas distribution are caused by supernovae. The 
holes are larger in Runs 3—6 since Eg^ is ten times higher in these runs compared to Runs 1 and 2. Consequently the gas discs of Runs 
1 and 2 are razor thin compared to Runs 3—6 which have a larger vertical gas velocity dispersion. 
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Run 3; t=0.32 Gyr 
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Figure 2. Star formation rate (SFR) history of the simulations. 
The solid, dashed, dot-dashed, dotted, triple dot-dashed and long- 
dashed lines represent Runs 1, 2, 3, 4, 5 and 6 respectively. 




10^^M(7) and the c oncentration parameter is set at c = 10 
iRobertson fc Kravtsov 2008. ). 

The stellar disc is assumed to follow an exponential sur- 
face density profile: 

Pd^ = ~A g-sech ( — )exp-( — ), (5) 

-i-iTZdUd Zd Hd 

where Md is the disc mass, Rd is the scale length and Zd is 
the scale height. Following suggestions from observations of 
M33 and carrying out several test runs at lower resolution, 
we decided to use Md = 4.0 x 10^ Mq, Rd = 1.4 kpc and 
Zd = 0.3 kpc. A summary of the numerical values used for 
our model parameters in the initial model setup is given in 
Table 1. We used 4 x 10^ stellar particles in our simulations 
giving each star particle a mass of 1 X 10* Mq. 

The gaseous disc is set u p following the method de- 
scribed in S pringel et al.l (|2005l ). The radial surface density 
profile is assumed to follow an exponential law like the stel- 
lar disc with a scale length, Rd,gas almost three times larger 
than the stellar disc. The initial vertical distribution of the 
gas is iteratively calculated to reach hydrostatic equilibrium 
assuming the equation of state calculated from our assumed 
cooling and heating function. The total gas mass is 2.0 x 10^ 
Mq and is equal to half the total stellar mass. Since we use 
2 X 10^ gas particles in each run, each SPH particle also has 
a mass of 1 x 10'' Mq. Note that this means both the star 
and gas particles have the same mass resolution in our sim- 
ulations. This is not by chance but required from our new 
modelling of star formation, feedback and metal diffusion 
shown in Section 2.1. We apply the s pline softening and 
variable softening length suggested by iPrice fc MonaghanI 
l|200iD for SPH particles. We set the minimum softening 
length at 198 pc, which corresponds to the required soft- 
ening for gas with solar metallicity and density of nn = 1 
cm~'^, i.e. the assumed density threshold for star formation. 
For the star particles, we applied a fixed spline softening 
of 198 pc softening length. We summarise our simulation 
properties in Table 1. Column 1 represents the NFW halo 
concentration parameter. The second column represents the 
virial mass; the third column is the virial radius; Columns 4 
and 5 represent the mass of the stellar and gas disc respec- 



Figure 3. Face-on (upper) and edge-on (lower) snapshot of the 
galaxy at t = 0.32 Gyr for stars (left column) and gas (right col- 
umn) colour coded by density for Run 3. The brightest regions 
are the densest. Large holes or bubbles are visible in the gas dis- 
tribution due to powerful SNe explosions. The large scale vertical 
turbulence this causes is visible in the edge-on gas distribution. 

lively. Columns 6 and 7 represent the scale lengths of the 
stellar and gas discs respectively and column 8 is the scale 
height for the stellar disc. 



3 RESULTS 

Using the galaxy simulation parameters shown in Table 1, 
we carried out six primary simulation runs, named Runs 1-6 
where we consecutively changed one extra model parameter 
between Runs 1-4. Run 5 is for comparison with Runs 1 and 

4 and Run 6 is our fine tuned 'best model'. Our varying 
model parameters include C, , Esn and Esw- Run 1 has low 
values for all three of our model parameters. We carried out 
a simulation run (Run 2) with C*=0.2, ten times higher star 
formation efficiency. We shall refer to this model as the 'high- 
C* model'. Another run (Run 3) was done with C,=0.2 and 
i?SN=1.0 X lO^'^ ergs. We call this run 'high-i?sN model'. In 
our fourth model run, we also increase Esw by a factor of 
ten. We call this model our 'highest-feedback model'. Run 

5 has the same high feedback parameters as in Run 4, but 
with a low C, value as in Run 1. Run 6 is identical to Run 5, 
except that we increase the value of C* from 0.02 to 0.05. We 
chose this value based on the results of the previous runs, 
believing this to be our best model. Table 2 summarises the 
simulation runs and model parameters. 

Fig. [T] shows face-on and edge-on snapshots of the stars 
and gas after 1 Gyr of evolution for our six runs. Large holes 
or 'bubbles' of varying sizes are clearly visible in the face- 
on gas distributions. These bubbles are produced mainly 
by SNe II explosions. The bubbles are much larger in Runs 
3—6 since the energy released per supernova explosion is ten 
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Figure 4. Three-dimensional visualisation of the early chaotic 
scenes in one of our high feedback simulations (Run 4). This snap- 
shot was taken with the observer embedded in the plane of the 
disc looking toward the centre of the galaxy. The yellow spheres 
represent the stars and the blue (diffuse) material is the gas. The 
colour coding is by density, so whiter or bluer regions represent 
more dense conglomerations of stars and gas respectively. The 
clumpy and highly inhomogeneous nature of the gas is clearly 
visible, as are the holes punched out by SNe. The larger stars are 
closer to the observer. 



Table 2. Model parameters 


Model 


C* 


-EsN 








(erg) 


(erg s-i) 


Run 1 


0.02 


1.0 X 10™ 


1.0 X lO^e 


Run 2 


0.2 


1.0 X 10™ 


1.0 X lO^'^ 


Run 3 


0.2 


1.0 X 10^1 


1.0 X 10^6 


Run 4 


0.2 


1.0 X 10^1 


1.0 X 10^'' 


Run 5 


0.02 


1.0 X 10^1 


1.0 X 10^7 


Run 6 


0.05 


1.0 X 10^^ 


1.0 X lO^'' 



times greater in these models. The large bubbles in these 
models may survive for longer than one galactic rotation. 
The models with high feedback and star formation efficiency 
parameters (Runs 3 and 4) produce the most turbulent ver- 
tical gas distribution as can be seen in the edge-on images. 
Conversely, the vertical gas distribution of the discs of Runs 

1 and 2 are razor thin in comparison. Clearly then SNe play 
a critical role in shaping the ISM properties. We shall look 
into this more quantitatively later. The vertical gas distribu- 
tion of Runs 5 and 6 are slightly less turbulent (compared to 
Runs 3 and 4) since these runs employ a significantly lower 
C, value. 

Fig. [2] shows the total SFR as a function of time for our 
different simulation models. From Fig. [5] we see that model 
Runs 1 and 5 have the lowest SFR. This is due to their low 
C« value. Runs 2, 3 and 4 have a global SFR varying between 
0.2 - 0.4 M0yr~^ over 1 Gyr of evolution. There are impor- 
tant differences between these models however. Model Run 

2 starts with a initially high SFR at 0.4 M0yr~^ and de- 
creases fairly consistently towards 0.2 M0yr~^ after 1 Gyr. 
This is because the high star formation efficiency of this 
model means an initially high SFR, but as time passes and 




5 L . . I . . . I . . . I . . . I . . J 

2 4 6 8 10 

Radius (kpc) 

Figure 5. Vertical gas velocity dispersion against radius seen at 
two different times in the same simulation (Run 3). The solid red 
line shows the result at an early and chaotic time whereas the 
blue dot-dashed line shows the velocity dispersion at much later 
times. During the large bubble phase, the velocity dispersion is 
unrealistically high. 

the gas reservoir from which the stars are formed quickly 
gets used up, the SFR correspondingly decreases. Runs 3 
and 4 initially have a lower SFR compared to Run 2. This is 
due to large outflows or 'bubbles' which drive out the gas at 
very early times. After this initial violent period the SFR in- 
creases and peaks after 0.5 Gyrs at 0.4 MQyr"^ after which 
it decreases slightly to 0.3 M0yr~^ by 1 Gyr of evolution. 
This means the higher energy produced from SNe and stel- 
lar winds initially slightly suppress star formation at very 
early times but that more gas becomes available for later 
star formation. The early time suppression of star forma- 
tion is slightly stronger in Run 4 compared to Run 3 since 
in this model we implement higher energy release from stel- 
lar winds as well as from SNe. Run 5 has identical feedback 
parameters to Run 4, with the only difference being that C, 
is low as in Run 1. Run 5 has a similar SFR history to Run 1, 
showing that the SFR is extremely sensitive to the star for- 
mation efficiency (C) as expected. Run 6 has roughly dou- 
ble the SFR of Run 5 as expected from its higher C* value. 
After approximately 0.5 Gyrs the SFR does not evolve sig- 
nificantly in all our simulation runs. Before 0.5 Gyrs the ISM 
can be extremely turbulent, especially for the high-feedback 
models, as can be seen in Fig. [3] which shows the 2-D face- 
on and edge-on stellar and gas distribution at an early time 
(0.32 Gyrs) for our high-iSsN model (Run 3). At this time, 
many large gas outflows are seen to occur. Fig. |4] shows a 
3-D visualisation of both the stellar and gas distributions 
of the galaxy at early times for Run 4. The clumpy and 
chaotic nature of the ISM is clearly apparent at these early 
times, with large gas outflows and bubbles visible. Fig. [5] 
shows that during the time when a large gas outflow occurs, 
the velocity dispersion of the gas component (measured as 
described below) is correspondingly very high. At a much 
later time when the disc is more settled (e.g. 1 Gyr), the ve- 
locity dispersion is correspondingly lower. Therefore 1 Gyr 
was chosen as a sufficiently long enough time to wait for the 
system to settle before carrying out our main analysis and 
comparison to the observations. 
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-10 12-1012 

Log Ih (^o/pc') 

Figure 6. Surface SFR density as a function of surface gas density for our six models at 333 pc resolution (grey plus signs). The surface 
SFR densities are calculated using the surface densities of young stars less than 0.5 Gyrs old. The diagonal dotted lines represent lines of 
const ant SFE, indicatin g the level of SgpR needed to consume 1, 10 and 100 per cent of the gas reservoir in 10* years, exactly as in fig. 
13 of iBigiel et al.l l l2008i ). Therefore the lines also correspond to constant gas depletion times of 10^'', 10^ and 10* yr r espectively. The 
red solid and dashed lines represent the mean and 1 sigma upper and lower bounds respectively for the lBigiel et al.l | |200S ) data between 
log Ejj values of and 1. The blue dashed horizontal and vertical lines in Run 2 are drawn to emphasize the region of low Eh and SsFR 
where there are barely any pixels present. The same region is marked in Run 4 to show that this model can successfully produce regions 
with both low Sh and SgpR- The location of these pixels are shown in Fig. [T] In this figure Ejj = Shi + Sh2- We also overplot the 
results obtained using azimuthally averaged ring annuli of width 1 kpc (black crosses connected by solid lines) . 
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Run 4 




Figure 7. Face on view of the stellar (left) and gas (right) distri- 
bution at t = 1 Gyr for model Run 4. The white stars represent 
the spatial locations of the pixels selected in Run 4 of Fig. |6] which 
contain both low surface gas and SFR densities (enclosed by the 
dashed blue lines and the axes). We can clearly see that all of 
these pixels are located in bubble regions where the gas density 
is low. 

Fig. [S] shows the relation between gas and SFR sur- 
face densities for our six sample runs where we have plotted 
the units to be ide ntical to the obser v ed sample of spira l 
galaxies in fig. 13 of iBigiel et al] (I2008D . iBigiel et all (|2008| ) 
show the surface SFR density against total hydrogen sur- 
face density. Taking advantage of our self-consistent chemo- 
dynamical modelling, we analyse the hydrogen fraction for 
each gas particle and self-consistently show the total hy- 
drogen surface density in Fig. [6] where in our notation: 
Eh = Shi + ■ The results were obtained using a pixel- 
by-pixel methodology using a square-cell resolution of 333 
pc X 333 pc and a smoothing length of two pixel box lengths. 
The surface SFRs are calculated using the surface densities 
of young stars less than 0.5 Gyrs old. We only include the 
pixels which have at least 10 young star particles contribut- 
ing to them. In Fig. [51 we also plot diagonal dotted lines 
of constant star formation efficiency (SFE), exactly as in 
IBigiel et al. I I2OO81' ) in order to make the comparison easier, 
indicating the level of Esfr needed to consume 1, 10 and 100 
per cent of the gas reservoir in 10* years. We also overplot a 
bes t fit line (sol i d red line) to represent the mean values for 
the IBigiel et a l. (2008) data between log Eh values of and 
1 (see e.q. 2 of .Bigicl et al. 2008 ) and show the 1 sigma upper 
and lower bounds (dashed red lines) . We can clearly see that 
the EspR for model Run 1 is too low at all values of Eh. This 
is because this simulation had a very low star formation effi- 
ciency. Increasing the SFE parameter (Run 2) does increase 
Esfr but at the expense of accurately following the slope 
of the observations. Our high-C* model (Run 2) therefore 
still fails to accurately follow the observed trends. The high- 
-EsN model (Run 3) and our highest feedback model (Run 4) 
better follow the observed slope and values for the SK law. 
Interestingly, the increased energy released from SNe explo- 
sions and stellar winds in these models help to increase the 
star formation rate at very low gas densities. The major dif- 
ference between Run 2 and Run 4 comes from the lack of 
pixels with low values for both Eh and Esfr in Run 2, un- 
like the higher feedback models such as Run 4. By selecting 
the pixels that fall into the bottom-left region (enclosed by 
the dashed blue lines) for Run 4 in Fig. [S] we investigated 
where these pixels were located and found that they predom- 
inantly lie in large holes blasted out by SNe (Fig. [7|. This 



suggests that SNe may play a vital role in shaping observed 
SK-laws, especially at low gas densities. Run 5 also employs 
high feedback parameters, but with a low C, value of 0.02 as 
in Run 1. Correspondingly, we see that Esfr is too low like 
in Run 1. In Run 6 we have increased the value of C* to 0.05. 
This represents our best compromised model. We see that 
increasing C* results in higher Esfr values and a better fit 
to the observations. We note that due to the relatively low 
mass of our simulated galaxy, which does not have a very 
high SFR or contain very dense regions, we cannot inspect a 
wide range of Eh parameter space and are therefore limited 
to making our comparisons in a narrow range of Eh . 

In Fig. [6] we also overplot the results obtained using a 
radial ring analysis method (black crosses connected by solid 
lines) whereby we have calculated the values of Eh and Esfr 
in azimuthally averaged ring annuli of width 1 kpc. We see 
that both the radial ring and pixel-by-pixel methods produce 
broadly consistent results, with the possible exception being 
that the radial ring method does not probe well the low Eh 
regime. Some of the data points in this regime represent ring 
annuli regions lying just outside the main star forming disc 
of our galaxy. 

In Fig. [8] we compare the kinetic energy density of the 
HI gas, Ek = (3/2)EHicrHi, against the SFR surface density 
as in iTamburro et al.l (|2009l ) . Here cthi = \/ -\- al where 
(Tth and (Jt are the velocity dispersions due to thermal and 
turbulent motions respectively. In addition, to mimic HI ob- 
servations, we only take into account the gas particles whose 
density is in the region 0.1 < pni < 50 cm~^ and have tem- 
perature less than lO** K. We also use the same units for 
the axes and show with solid lines the supernova energy 
input for different SNc effici ency v a lues, e sN = 1, 0.1 and 
0.01 exactly as in Tamburr o~et all l|2009l) . Their observa- 
tions show a linear correlation between and Esfr, with 
typical esN values between 0.1 — 1. At higher Esfr regimes, 
the -Efc slightly falls off. Comparing against our simulated 
galaxies with their Fig. 4, we see that Runs 3—6 best follow 
the observations. Runs 1 and 2 fail to capture the observed 
trends, especially, i?fe is lower at —9.5 < Log Esfr < —8.5 
compared to Fig. 4 of ITamburro et al] (|2009| ). 

We also analysed the azimuthally averaged radial pro- 
file for the HI velocity dispersion, by calculating the av- 
erage value of (thi (following the same method as above) 
within ring annuli of radial width 1 kpc. Observations of 
spiral galaxies show a radial decline for the HI velocity dis- 
persion with a value of ~ 20 km s ~^ in the inner region , 
declining to ~ 10 km s"'^ by r25 l|Tamburro et al.ll2009l ). 
Fig. |5] shows that our higher feedback models (Runs 3 and 
4) closely follow the aforementioned observed trends. All of 
our models had similar cthi values of aro und 10 km sT^ b y 
r25 similar to what was found in iTamburro et al.l (|2009l ). 
With the exception of Run 1, all models showed a radial 
decline of cthi. The decline was small for our high-C* model 
(Run 2) and more noticeable for our higher feedback mod- 
els (Runs 3 and 4). This is understandable since the higher 
feedback models produce larger gas outflows and bubbles 
which increases cthi. Interestingly, our model with high i5sN 
and -Esw (Run 4) did not have a noticeable difference with 
our high-i7sN model (Run 3), proving that it is the super- 
nova energy which mainly drives the turbulence in the ISM 
in our models. Further evidence of this lies in the compari- 
son of the velocity dispersion of Runs 1 and 5, with Run 5 
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Figure 8. Pixel-by-pixel scatter plot of kinetic energy vs. surface star formation rate for our six model galaxy runs (grey plus signs). 
The solid lines of unity slope represent the supernova energy input for different supernova efficiency values, csn = 1,0.1,0.01 (top to 
bottom). 
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Figure 9. Gas velocity dispersion versus radius for our six models. The 'higher-feedback' with high-C* models (Runs 3 and 4) have 
larger velocity dispersion. Our results show that the velocity dispersion is primarily sensitive to our SNe energy (Ssn) parameter. 



having a noticeably larger (thi- Despite this, cthi for Run 5 is 
still to o low compared with observations of iTamburro et al] 
Therefore in Run 6, we slightly increased the value 
of C* and found that cthi also increased to a level which is 
consistent with the observations. 



4 SUMMARY AND CONCLUSIONS 

We aim to build the most accurate and self-consistent nu- 
merical model for late-type spiral galaxies to date. To this 
end we have been developing our numerical simulation code 
GCD+. We include new star formation and feedback recipes 
which allow our simulations to follow the effects of powerful 
SNe explosions on galaxy evolution. 

In this study we set up an isolated M33-sized disc galaxy 
consisting of gas and stellar discs with no bulge compo- 
nent, in a static dark matter halo gravitational potential. 
We carried out six main simulation runs in which we varied 
three of our model parameters; the star formation efficiency 
(C), the thermal energy released per SNe explosion (Esn) 
and the thermal energy released per unit time from stellar 
winds (iJsw)- We studied the effect of varying these parame- 



ters on the Schmidt-Kennicutt Law and vertical gas velocity 
disp ersion versus radius rela tion comparing to obse rvations 
from lsigiel et all (|2008l ) and lTamburro et al.1 (|2009l ) respec- 
tively. We find that our models with higher energy feedback 
(Runs 3—6) more closely resemble the observations. Power- 
ful supernova explosions cause the formation of large holes 
or bubbles in the gas distribution in the galaxy and stir the 
ISM turbulence and are thus essential for producing the ob- 
served gas velocity dispersions. The gas velocity dispersion 
is also affected (albeit less strongly) by the star formation 
efficiency. Our model with high energy feedback from SNe 
and stellar winds and intermediate star formation efficiency 
(Run 6) most closely resembled the observations. This run 
represented our 'best compromised model' since we had fine- 
tuned the star formation efficiency and feedback parameters. 
In our pixel-by-pixel analysis of the star formation law, we 
found that our model with high star formation efficiency and 
low energy feedback (Run 2) produced a KS relation with an 
unrealistically steep slope. Models with low energy feedback 
had a distinct lack of pixels with both low Eh and Esfr, re- 
sulting in a failure to accurately reproduce the observed star 
formation law. By increasing the energy feedback parame- 
ters, i5sN and -Esw (as in Run 4) we managed to rectify this 
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problem and found that the pixels which had low values for 
both Eh and Esfr were predominantly located inside bub- 
ble regions, which were more numerous and larger in size for 
the higher feedback runs. Our results therefore suggest that 
SNe can play a pivotal role in shaping the star formation 
law. We used both a pixel-by-pixel and radial ring method 
to analyse the star formation law in our galaxies and found 
that overall, both methods gave broadly consistent results. 

In this paper we employed static potentials for the dark 
matter halo and stellar disc. The main benefits of not using 
a live halo include the prevention of numerical efi'ects such 
as artificial disc heating due to the scattering of baryonic 
particles by more massive dark matter particles which may 
be an issue in lower resolution cosmological simulations. Our 
numerical simulations thus do not sufi^er from these global 
disc instabilities. However, we und erestimate the pert urba- 
tion from accreting s atellites (e.g. Purcell et al.ll201ll ) and 
globular clusters (e.g. IVande Putte et al.ll2009l '). 

In this paper we demonstrate that reproducing both the 
SK relation and the observed gas velocity dispersions puts 
stronger constraints on the sub-grid modelling for galaxy 
evolution simulations. Encouraged by the success of this 
work, we are currently in the process of running higher 
resolution simulations. This should result in us being able 
to probe deeper into the origins of the Egas — Esfr rela- 
tion. It should also mean that our models will be less sensi- 
tive to the employed value of the star formation efficiency, 
;iven we could res olve molecular cloud formation processes 
Saitoh et al I2OO8'). After calibrating our code at this higher 
resolution, we then aim to run new extremely high resolution 
cosmological simulations. The fully self-consistent nature of 
the new higher resolution cosmological simulation will be a 
powerful tool to unravel the mysteries behind the formation 
and evolution of late type spiral galaxies like M33. 
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